﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using NUnit.Framework;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics;


namespace QuantitativeIndicator.Tools
{
    class MathLib
    {
        //求最小二乘法线性回归的系数
        public double[] regress(double[] Y, double[] X)
        {
            double[,] Y1 = new double[Y.Length, 1];
            for (int i = 0; i < Y.Length; i++)
            {
                Y1[i, 0] = Y[i];
            }

            double[,] X1 = new double[X.Length, 2];
            for (int i = 0; i < X.Length; i++)
            {
                X1[i, 0] = 1;
                X1[i, 1] = X[i];
            }

            Matrix mY = Matrix.Create(Y1);
            Matrix mX = Matrix.Create(X1);
            Matrix mResult = mX.Solve(mY); ;
            return mResult.GetColumnVector(0);
        }

        //求得线性回归残差
        public double[] getResidualError(double[] Y,double [] X)
        {
            double[,] Y1= new double[Y.Length,1];
            for (int i = 0; i < Y.Length; i++)
            {
                Y1[i, 0] = Y[i];
            }

            double[,] X1 = new double[X.Length, 2];
            for (int i = 0; i < X.Length; i++)
            {
                X1[i, 0] = 1;
                X1[i, 1] =X[i];
            }

            Matrix mY = Matrix.Create(Y1);
            Matrix mX = Matrix.Create(X1);
            Matrix mResult = mX.Solve(mY);
            Matrix mResidualError = mY - mX.Multiply(mResult);

            return mResidualError.GetColumnVector(0);
        }
        //求残差的平方和
        public double getResidualErrorSquare(double[] Y, double[] X)
        {
            double [] mResidualError=this.getResidualError(Y,X);
            double sumR = 0;
            for (int i = 0; i <mResidualError.Length; i++)
            {
               sumR=sumR+ Math.Pow(mResidualError[i], 2);
            }
            return sumR;
        }
        /// <summary>
        ///    求得一元线性回归的先关系数
        /// </summary>
        /// <param name="Y">因变量</param>
        /// <param name="X">自变量</param>
        /// <returns></returns>
        public double getCorrelationCoefficient(double[] X, double[] Y)
        {
            double r = 0;
            r = (X.Length * this.multiplierSum(X, Y) - this.sum(X) * this.sum(Y)) / Math.Sqrt(X.Length*this.squareSum(X)-Math.Pow(this.sum(X),2))/Math.Sqrt(X.Length*this.squareSum(Y)-Math.Pow(this.sum(Y),2));

            return r;
        }
        
        /// <summary>
        /// 一维数组和
        /// </summary>
        /// <param name="X">一元数组X</param>
        /// <returns></returns>
        public double sum(double [] X)
        {
            double sum = 0;
            for (int i = 0; i < X.Length; i++)
            {
                sum = sum + X[i];
            }
            return sum;
        }
        /// <summary>
        /// 一维数组平方和
        /// </summary>
        /// <param name="X"></param>
        /// <returns></returns>
        public double squareSum(double[] X)
        {
            double sum = 0;
            for (int i = 0; i < X.Length; i++)
            {
                sum = sum + X[i] * X[i];
            }
            return sum;
        }
        /// <summary>
        /// 两数组相乘和
        /// </summary>
        /// <param name="X"></param>
        /// <param name="Y"></param>
        /// <returns></returns>
        public double multiplierSum(double[] X, double[] Y)
        {
            double sum = 0;
            for (int i = 0; i < X.Length; i++)
            {
                sum = sum + X[i] * Y[i];
            }
            return sum;
        }
        public static void test()
        {
            MathLib mathLib = new MathLib();
            double[] y = { 1, 3.4, 5.6 };
            double[] x = { 2, 4.5, 5.4 };
            double[] t = mathLib.getResidualError(y, x);
            double sR = mathLib.getResidualErrorSquare(y, x);
            double r=   mathLib.getCorrelationCoefficient(x,y);
            double d = r * r;
            Console.WriteLine("");
        }
    }
}
